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In the present paper we give a brief summary of some recent theoretical advances in the treatment 
of inhomogeneous fluids and methods which have applications in the study of dynamical properties 
of liquids in situations of extreme confinement, such as nanopores, nanodevices, etc. The approach 
obtained by combining kinetic and density functional methods is microscopic, fully self-consistent 
and allows to determine both configurational and flow properties of dense fluids. The theory predicts 
the correct hydrodynamic behavior and provides a practical and numerical tool to determine how the 
transport properties are modified when the length scales of the confining channels are comparable 
with the size of the molecules. The applications range from the dynamics of simple fluids under 
confinement, to that of neutral binary mixtures and electrolytes where the theory in the limit of 
slow gradients reproduces the known phenomenological equations such as the Planck-Nernst-Poisson 
and the Smolochowski equations. The approach here illustrated allows for fast numerical solution 
of the evolution equations for the one-particle phase-space distributions by means of the weighted 
density lattice Boltzmann method and is particularly useful when one considers flows in complex 
geometries. 

PACS numbers: 


I. INTRODUCTION 


One of the goals of statistical physics is to bridge different levels of descriptions of matter accounting for macroscopic, 
mesoscopic and microscopic levels. A complete description at the finest scale makes in principle possible to derive the 
g^erties at the grosser scales using a reduction process usually carried out only through a series of approximations 

In Thermodynamics and Hydrodynamics, which both refer to the macroscopic level, a different strategy is adopted 
and instead of using a reduction process, one deduces relations between a reduced number of macroscopic variables 
such as tenxperature, density, energy, fluid velocity, etc. from experimental observations and/or phenomenological 
arguments [j, 0 . 

In principle one may obtain information about large assemblies of interacting molecules via molecular dynamics 
(MD) numerical simulations, but the method requires considerable computational time and computer memory es¬ 
pecially in systems where one is interested to explore properties relative to species having very low concentration. 
Alternatively, the theoretical methods which can be used to study the dynamical behavior of fluid solutions are the 
DDFT (dynamic density functional theory), which describes an overdamped dynamics typical of colloidal behavior 
0-0 the WDLBM (weighted density Lattice Boltzmann method) which describes structural and thermodynamical 
features together with the inertial dynamics of liquids 0, [l0| . 

DDFT derives a closed evolution equation for the particle density, n, by considering the local conservation law 
for the number of particles and relating the current density to the gradient of the density itself by using the rules 
governing the microscopic dynamics of the system or some phenomenological argument. The existence of a simple 
relation, useful in practical applications, between current and density is not a priori guaranteed and in some cases, 
depending on the nature of the microscopic dynamics, it might be more convenient to explicitly consider, in addition 
to the equation for n, the balance equations for the momentum density and energy density variables. The use of this 
larger set of variables to describe the state of the system is characteristic of hydrodynamics or of the kinetic approach, 
based on the knowledge of the phase space distribution function. 

As far as hydrodynamics is concerned the application of the Navier-Stokes equation at the nanoscale has been 
questioned, because the presence of surfaces or the confinement of the particles in very tiny regions may affect 
drastically the local relations between hydrodynamic variables. Kinetic theory is the natural tool to treat dynamical 
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effects at the microscopic level, but requires a large effort because one needs to consider the phase space (r, v) instead 
of the standard coordinate space. However, recent studies have shown that it is possible to conjugate the features of 
the DFT with those of the discrete numerical technique introduced about two decades ago to solve kinetic equations, 
namely the Lattice Boltzmann method (LBM). The standard LBM suffers from a lack of physical realism in the 
treatment of the interactions and of thermodynamic consistency. It is not fully microscopically motivated, but rather 
is inspired by a top-down approach. It describes hydrodynamic effects, but has a trivial or ”ad hoc” thermodynamics 
corresponding to structureless fluids. Shan and Chen [ij and He et al. [l^ proposed independently modifications of 
the LBM aimed to include a better thermodynamic treatment of the fluids, however none of these approaches were 
able to capture the structural features of fluids which are among the merits of the DDFT theory. 

The present paper is organized as follows in sect. |n]we give a short account of the derivation of the kinetic equation 
for the single particle phase space distribution function, in section Hill we consider the case of a colloidal solution where 
there is a scale separation between the colloidal particles and the solvent particles and is possible to represent the 
system as an assembly of Brownian particles governed by overdamped dynamics. In section HVl we consider instead 
the case where solute and solvent particles have similar physical properties and the solvent has to be treated explicitly. 
The resulting dynamics is inertial and one has to take into account adequately the hydrodynamic modes. In section 
m we give explicitly the extension of the kinetic equation to the study of neutral and charged mixtures, while in 
section El we briefly illustrate how the kinetic equation can be solved numerically using tools borrowed form the 
lattice Boltzmann method. Finally in section rVHI we present a brief summary. 


II. KINETIC DESCRIPTION OF INHOMOGENEOUS FLUIDS 

Let us consider an assembly of N particles, of mass m and positions and velocities v^, mutually interacting through 
pair potential U{ri — rj) and subject to external forces Fexti'Ti)- The system evolves according to the following set of 
equations: 


drj 

dt 


= Vi 

= -Vr.H(u)- ^ Vr 


U{\r^ - Vjl) - mjv, + ^^{t) 


( 1 ) 


where the last two terms in o represent the coupling to a stochastic heat bath at temperature T, characterized 
by a stochastic white noise forcing whose noise amplitude is related to the friction constant, 7 , by the Einstein 
fluctuation-dissipation relation (s)) = 2^mkBTSij6°‘^6{t — s), where ks is the Boltzmann constant. The noise 

represents the effect of microscopic degrees of freedom not explicitly accounted for in the interaction terms. 

For our purposes it is more convenient to switch from the coupled set of stochastic differential equations (P) to 
the description based on the N-particle probability density distribution /Ar({r, v}, t), where {r,v} indicates a 67V 
dimensional phase space point, which evolves according to the Kramers-Fokker-Planck 0 : 


^/v({r, v}, t)+Y^ • Vr, - ( 


Vr.E(u) ^ Vr,[/(|x, -x,|))-^j/jv({r,v},t) 


m 




= 7 




keT 
m - 


/v({r,v},t) 


( 2 ) 


The probabilistic description represented by eq. (P is the result of an ensemble averaging of the trajectories over 
a noise ensemble and over initial conditions in the case of damped stochastic dynamics (7 > 0 ) typical of a colloidal 
suspension. In the case Hamiltonian dynamics characteristic of an atomic liquid, where 7 = 0, eq. (P reduces to the 
Liouville equation. 

The information contained in /^v is fully microscopic and describes the probability density of the microstates of 
the system. However, the distribution /^r is a very complicated object to handle and provides a representation 
redundant for practical purposes. One then seeks a contracted description from the 67V dimensional F-phase space 
to a 6 /r-dimensional space, that is, from the phase space distribution of TV particles to the one particle phase space 
distribution, /(r, v,7). In order to contract the description from the F to the ^-space, one considers the marginalized 
s-particle distribution functions, fs, constructed by integrating over the 2d degrees of freedom of (TV — s) particles 


TV! 

{N-s)l 


N 

dVndVn /Ar(ri, ..,rs,rs+i..,rAr, Vi,.., v,,, v^+i.., vat, T). 

n=s+l 


/,,(ri,..rs,vi,..,v,,,t) 


( 3 ) 
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It is straightforward to derive a set of of coupled equations connecting the distribution functions of different orders s 
among themselves, the so called Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy 1^. The one particle 
distribution is / = /i, which is the object of the kinetic theory, is a useful bridge between the many-body molecular 
dynamics models and the continuum mechanics models. The evolution equation for / , the so called kinetic equation, 
contains a non linear term, named collision term, which takes into account the effect of interactions in an approximate 
fashion. 

At the bottom end of the hierarchy, that is integrating equation ((H])) over {N — 1) particle’s coordinates and 
velocities, one obtains the exact evolution equation for the one particle distribution: 


Equation (l4|) contains the streaming terms in the l.h.s., and the interparticle interaction term, 17, the r.h.s. and a 
coupling to the stochastic heat bath : 


= 7 

For continuous pair potentials, U, one can write: 

1 d 


m dv'^ dv 




• J dr'y dv72(r,v,r',v',7)Vrf/(|r- r'l). 


(5) 

( 6 ) 


Instead, in the case where U contains an hard-sphere contribution, VrG is singular and one has to treat 17 differently: 
the non singular long range part of the potential can still be handled as in eq.®, while the singular piece has to be 
dealt by special methods of kinetic theory [Il[i3. In fact, the dynamics of hard spheres is no longer described in 
terms of forces since the trajectories of the particles abruptly change when a pair comes into contact and exchanges 
momentum according to the law of elastic collisions. By using the so-called pseudo-Liouville operator approach one 
reformulates the dynamics in terms of binary collision operators rather than forces. This method leads to the revised 
Boltzmann-Enskqg kinetic equation describing a system of dense hard-spheres and was introduced rigorously by Ernst 
and van Bejeren jT3, S] about 40 years ago. In the following we shall use this approach to treat the collision term, 
whose explicit form is given by eq. (USD, when hard sphere interactions are involved. 

In a nutshell, the underlying assumption is that the N-particle distribution function /at at all times takes exactly 
into account the non overlap of any pair of spheres. Instead, the velocity dependence of fjy is factorized into a product 
of single-particle velocity distribution functions and the resulting kinetic equation for the one-particle distribution 
function, /(r, v,7), involves the two-particle distribution function, 7 . As suggested by Enskog, /2 is expressed as the 
product of one particle distribution functions times the positional pair correlation function g{r,r\t): 


/ 2 (r, V, r', v', t) Ri /(r, v, t)/(r', v', t)g{r, r', t) 


The approximation employing the product of two one-particle distribution functions is called ’’molecular chaos” 
hypothesis and it disregards the correlations of the velocities of the two colliding particles, prior to the collision. 
Importantly, it decouples the evolution of /(r, v,7) from the evolution of the higher order multiparticle distribu¬ 
tion functions. Many-particle correlations are however retained through the structural information contained in the 
positional pair correlation function function ^(r, r'). As an approximation, we shall assume that g{r,r') is a non 
local function of the profile n{r,t), depends on time only through the density profile and has the same form as in a 
nonuniform equilibrium state whose density is n{r,t). 

In order to describe the properties of molecular fluids we need the hydrodynamic framework to consider five field 
variables as the minimal ingredients to describe the local (macro)state of a simple fluid. These fields are the first 
velocity moments of /(r, v, t) and correspond to the number density, local velocity and local temperature of the fluid: 


/ n{r,t) 
n(r,7)u(r,t) 

\ lkBn{r,t)T{r,t) 





(7) 


The hydrodynamic fields obey the following balance equations, obtained by multiplying the kinetic equation by the 
first Hermite tensorial polynomials, proportional to {l,v, (v — u(r, 7))^}, and integrating over the velocity : 


5tn(r,7)-I-V • (’^(i’,7)u(r,7)) = 0 (8) 

mn{r,t)[dtUj{r,t) + Ui{r,t)ViUj{r,t)] + ViPl;f\r,t) + n{r, t)V jV (r) - n(r, 7)$^ (r, 7) = bf\r,t) (9) 

^kBn{r,t)[dt + Ui(r, 7)Vi]r(r, 7) -b {Y,t)ViUj{v,t) -b ViqP{Y,t) - Q{r,t) = b^'^\r,t). 


( 10 ) 
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We have introduced the kinetic component of the pressure tensor and the heat flux vector 


= /dvf ) /(r,v,t). 


The ’source’ terms stem from the heat bath 




6(2) (r, t) j - y 1, f (v - u)2B(r, v, t) j “ 1, < > -3^r(r, t) 

Finally, one considers the terms due to irrterparticle interactions: 

n(r,t)$i(r,t) =m j dvU.{v,v,t)vi = ^{r, t) 


( 11 ) 


( 12 ) 


(13) 


(C) 

It can be shown that 4>i is related to a spatial derivative of the potential part of the stress tensor Pb while 

Q(i’:l) = y y'lvn(r,v,t)(v - u)2 (14) 

is a combination of the gradient of the collisional contribution to the heat flux vector and the pressure tensor times 
the strain rate: 


Q{r,t) = -Vjg|(r,t) -Pb(r,t)VjUi(r,t). 

Both <i>i and Q vanish in uniform systems. 

The difficulty with eqs. ([8l) - (ll0ll is not only due to the non linearity stemming from the interaction terms, but 
also from the projection of the kinetic equation onto the velocity space ( spanned by the various Hermite tensorial 
polynomials ). In fact, the procedure eliminates the velocity dependence, but on the other hands generates an infinite 
hierarchy of equations for the evolution of the velocity moments of /(r, v, t) due to the coupling determined by the 
presence of the convective term v • Vr/. Notice that and q^^^ cannot in general be expressed in terms of the 

five hydrodynamic moments. To obtain a closure, one must truncate this hierarchy at a given level, which requires 
the approximation of higher moments. For example, for the standard truncation at the velocity (n = 1) level (i.e. the 
same level of description as the Navier-Stokes equations), one must control terms (in appropriate units) of the form 
/ dv (v 0 V — /)/(r, v,t), where I is the 3x3 identity matrix. One way to deal with such a problem is to introduce 
phenomenologically motivated relations which allow to express higher order velocity moments in terms of the five 
basic moments, these are the so called constitutive equations for the momentum and heat flux. In the case of systems 
described by (under)damped dynamics it is possible to eliminate all moments higher than the zeroth moment, in favor 
of n(r, t) by multiple time scale expansion | 1 ^ . 

However, as recently proposed by our group it is possible to solve directly the kinetic equation even in the presence 
of non trivial collision terms by using the Weighted Density Lattice Boltzmann Method (WDLBM) which combines 
the structural features of the DDFT with the requirements of hydrodynamics and also gives transport coefficients in 
self-consistent fashion, see section Hvl 

To proceed further we must give explicitly the form of D and discuss how to treat it. In order to obtain a suitable 
description of the fluid behavior, reproducing at least the qualitative feature of the equation of state and the density 
dependence of the transport coeSicients the interaction must be separated into short range repulsion (assimilable to 
hard spheres) and long range attractive contributions: 

D(r, V, t) = ftrep{r, V, t) + Dottr(r, V, t) 

A good choice is to assume for flrep the Revised Enskog Theory (RET) form for hard spheres pill: 

Drep(r, V, t) = J dv2 J (ifi0(6 • Vi2)(6 • V12) x 

{52(1’, r - fi 6 )/(r, v')/(r - 66, V2) - g2(r, r + 66 )/(i', v)/(r 66, V2)} (15) 

where the primes indicate post-collisional velocities v' = v — (6 • Vi 2)6 and V 2 = V 2 -|- (6 • Vi 2 ) 6 . determined by the 
conservation of total momentum and total energy in an elastic collision; a is the hard-sphere diameter and 6 the unit 
vector directed from the center of sphere 1 to the center of sphere 2. Interactions are non-local and momentum and 
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energy can be transferred instantaneously across finite distances when the spheres collide An adiabatic approximation 
for 52 (r, r', t) is chosen by which 52 (r, r', t) is assumed to be given by the equilibrium pair correlation function of the 
system when its ensemble averaged density is given by n(r, t). In addition since the exact form of the inho mog eneous 
pair correlation 32 is not known, one uses the Fischer-Methfessel prescription to construct it locally [^. This 
ansatz is common both to the DDFT and the WDLBM methods. The calculation of the collision integral involves the 
inhomogeneous pair correlation function, 32 , at contact which is not known. One resorts to the Fischer and Methfessel 
prescriptionji^, which assumes that the functional form of the inhomogeneous (72 (r, r + crk) is obtained by replacing 
the density n(r) by the associated coarse grained density h(r). The latter is the average of n(r) over a sphere of 
diameter a centered at r: 


n(r) =[ dr'n{r + r') 9 {^ - \r - r'\) 

6 ^ 7 2 

The following rule gives the pair correlation function at contact: 

52(r,r + crk) =5^“"=(j7(r+ika)) ( 16 ) 

where the local average pac king fraction fj is fi{r) = ^n{r)a^ and the explicit expression of is provided by the 
Carnahan-Starling equation [ij : 


Finally, the attractive contribution to the collision term is written in the mean-field form 

f 2 attr(r,v,t) = ^Vv • y dr' J dv'f{r,v,t)f{r',V,t)WrU{\r-r'\) (18) 

where configurational correlation are disregarded having set g{r,r',t) = 1 in this formula. Due to such crude approx¬ 
imation the attractive term contributes to the equation of state, but does change the transport coefficients, with the 
notable exception of the diffusion coefficient. 


III. OVERDAMPED DYNAMICS 

In this section we shall specialize to the over damped dynamics and obtain the DDFT equation of evolution. Classical 
density functional theory (DFT) has been used with great success to investigate the structural and thermodynamic 
properties of inhomogeneous classical fluids [2l|. It represents a relevant generalization of the highly successful and 
rigorous method originally introduced about 50 years ago by Hohenberg and Kohn in quantum many body theory. 
The widespread use of the DFT approach is due to the fact that its fundamental entity, the Helmholtz free energy, 
A[n], is an intrinsic functional of the local density of molecules n(r’), i.e. is independent of the external potential 
to which the fluid is subjected, in other words A is a universal functional once the interactions among the particles 
and their properties are specified. The theory states that for a fixed external potential V{r), fixed temperature and 
chemical potential, there is a unique equilibrium density profile n(r), which can be found by minimizing the Grand 
potential functional f2[n] = F[n\ -|- f n(r)[H(r) — fj]dr with respect to n. Although F[n] is known exactly only in few 
particular cases, fairly good approximations have been devised, so that the method is versatile and generally applicable 
with success to study the properties of real systems under a variety of thermodynamic and geometric conditions. A 
variety of problems has been solved, ranging from adsorption, phase transitions at surfaces, wetting phenomena, fluids 
confined in narrow pores, theory of freezing, depletion forces, etc. 

Under appropriate conditions such as those realized in colloidal solutions, recent studies have shown that a dynamical 
extension of the DFT, the so called dynamical density functional theory (DDFT) may provide a valid description of 
the main features of the approach towards equilibrium. In order to establish the DDFT equation of evolution within 
the present framework one must consider the role of the two terms in eqs. (EJ-dlOl), accounting for the drag 

force proportional to the particle velocity and the random stochastic acceleration due to the solvent atoms impinging 
on the molecules. These terms determine a fast equilibration process of the momentum current nu and of the local 
temperature T towards their stationary values when the friction 7 is large (overdamped limit). In this case, one may 
neglect the inertial term dvj jdt in eq. ([T]) and derive the N-particle Smoluchowski equation associated with such first 
order system, a time dependent partial differential equation for the distribution function of the V particle positions. 
This is the starting point of the strategy followed by Archer and Evans Q to derive their version of the DDFT by 
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integrating out the coordinates of {N — 1) particles. The presence of the friction terms is of great help in reducing the 
infinite hierarchy, whose five equations appearing in eqs. diD-dini) are only the first members, to the single equation 
for the density appearing in the DDFT equation. 

The rigorous mathematical procedure to derive this result employs the multiple time scale analysis , which 

exploits the time scale separation between the zeroth mode associated with the density fluctuations and the remaining 
modes, which also include the momentum and energy fluctuations. Due to the action of the heat bath these are fast 
relaxing modes because the velocities of the particles rapidly relax towards the equilibrium distribution, in a time of 
order of the inverse friction time r = I/ 7 . Since the momentum and energy of the colloidal particles are not conserved, 
their currents become rapidly “slaved” to the density, i.e. the evolution is completely determined in terms of n(r, t) 
and its derivatives. The only relevant evolution on timescales larger than t = I /7 regards the spatial distribution of 
the particles. This is the reason why the DDFT gives a sufficiently accurate description of colloidal systems. One can 
determine the particle current, J(r,t) = n(r, t)u(r, t), by imposing the vanishing of the inertial terms in eq. ([H|) and 
using eq. 0 , so that from the momentum balance the following approximate equality holds: 

+ ViP'>f\Y,t) + n(r,t)VjD(r) « -m'yn{r,t)uj{r,t). (19) 

Finally, substituting nuj into the continuity equation ([5]) one finds: 


dt my 


+nir,t)V^V{r,t) 


( 20 ) 


Neglecting dissipative contributions (viscous and thermal conduction, see eq. (l26l) below) to the pressure tensor we 
can write pj;^\r,t) = kBTn(r,t)5ij or equivalently 




(i’,i) 


n(r,t)Vi 


ideal [^] 

6n{r, t) 


where 


Pideai[n] = ksT J drn{r,t)(\nn{r,t) - 1 ). 
is the ideal gas contribution to the free energy. 

Similarly one finds the following relation between the excess pressure and the excess free energy. Pint (see eq. (|27ll 
below and references): 




(C) 


(i’,i) 


n(r,t)V, 


int [^] 

6n{r, t) 


Collecting together one obtains the following equation 


9 n(r, t) 
dt 


—Vin(r,t) 
rri'y 


f ^Pideal \yi\ 

* V 5n{r,t) 


^Pint [^] 

6n{r, t) 



( 21 ) 


One can recognize that equation (EU for the density corresponds to the DDFT equation. The details of the 
derivation of the DDFT formula can be found in ref. Such a formula has a kinetic derivation, but one can 

make contact with the equilibrium density functional theory by considering the evolution of the statistical average 
of the microscopic instantaneous density over an ensemble of identical copies of the original system, each copy being 
characterized by a different realization of the noise. Such a procedure, besides being easily realizable in a numerical 
simulation, leads to a simplification of the equation (12111 and to the interpretation of the drift term as the functional 
gradient of the derivative of the free energy functional with respect to the density, that is the local chemical potential 
gradient. 

The crucial assumptions made in deriving eq. (1^ are the following: a) the density n(r, t) is the average of the 
microscopic instantaneous density over the realizations of the random noise and is therefore a smooth function of the 
coordinate r. b) The functional F[n] is a function solely of the density n{r,t). c) The instantaneous two-particle 
correlations are contained in Pint [n] and are approximated by those of an equilibrium system having the same density 
profile as the system at time t. This is the so called adiabatic approximation. Since in DDFT the driving force 
towards the steady state is given by a derivative of the free energy it cannot contain contributions from dissipative 
forces. Such an aspect is at variance with the kinetic approach where, as shown below, frictional and viscous forces 
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appear naturally in the stress tensor. From eq. (1211) one can derive an H-theorem stating that the free energy never 
increases during the relaxation process. 

Notice that the theory has selected a particular reference frame, where the solvent is at rest. Such a choice breaks 
the translational invariance of the system, which for this reason cannot support sound waves. 

The DDFT has been applied to several problems and extended to treat non simple fluids and more complicated 
interactions such as the Hydrodynamic interaction [ 2 ^. It can be extended to describe charged and uncharged fluid 
mixtures and applies whenever the system dynamics is a diffusion relaxation process. Therefore it can be used to 
describe the dynamics of ions in a aqueous solvent, if the structure of the latter is not important for the problem. 

The overdamped dynamics illustrated above is the extreme limit where ycr » vt, where vt is the thermal velocity. 
However, one can investigate corrections for finite values of 7 [ 2 ^,[ 2 ^. These corrections can be computed by an inverse 
friction expansion or by multiple time scale methods and give rise to a more complex DDFT equation as shown in 
refs. [ 2 ^, [23 , which takes into account momentum and energy fluctuations. 


IV. INERTIAL DYNAMICS 

The inverse friction expansion does not help in the case of molecular fluids where the inverse friction parameter 1 /y 
diverges, and only internal dissipation mechanisms are at work. We consider now the case with the 7 = 0 ( 6 ^*^ = 0 
). A salient feature of molecular liquids is their ability to support hydrodynamical modes, since particle number, 
momentum and energy are locally conserved. The set of evolution equations for the five hydrodynamic variables do 
not form a closed system, unless one assumes some phenomenological constitutive relations between the gradients of 
Pij and q (the sums of the kinetic and collisional parts) and the fluxes. 

Using the kinetic equation it is possible to approximate the non linear interaction term and obtain its explicit 
dependence on the small set of hydrodynamic modes and not on the full distribution function /(r, v,t) obtaining a 
fast numerical solution of the equation. Such a step is very important in our treatment because it allows to reduce 
enormously the computational effort and obtain a fast numerical solution of the kinetic equation. Such a reduction 
must preserve the translational invariance of the system, i.e. must not select the reference frame where the solvent is 
at rest 

A practical treatment of the RET collision operator was suggested by Dufty et al [13,mi, who proposed to separate 
the contributions to D stemming from the hydrodynamic modes from those from the non-hydrodynamic modes. Such 
a goal is achieved by projecting the collision term onto the hydrodynamic subspace spanned by the functions { 1 , v, v^} 
and onto the complementary kinetic subspace: 


with 


Phydro^ — 


kBT{r,t) 

where $ and Q are given by eqs. CSD and (na, respectively and 


^ — T^hydro^ “ 1 “ {J T^hydro)^ 


t) (v - u) • n(r, t)$(r, t) + - t) 


(j)M{r,v,t) = [ 


27rfcBT(r, t) 


] 3/2 


exp 


ikBT{r,t) 

m(v — u)^' 
2kBT{r,t) 


( 22 ) 

(23) 


is the local Maxwellian. As far as the projection of ft onto the non-hydrodynamics sub-space Dufty and coworkers 
approximated it by a phenomenological Bathnagar-Gross-Krook (BGK) prescription, which preserves the number 
of particles, the momentum and the kinetic energy, thus fulfilling the physical symmetries and conservation laws of 
the fluid: 


(7 - 'Phydro)ft ^ -V /(r, V, t) - n(r, t)(j)M (r, v, t) 


(24) 


This part of the collision operator contributes to determine the values of the transport coefficients which turn out 
to be functions of phenomenological parameter v, a phenomenological collision frequency chosen as to reproduce the 
kinetic contribution to the viscosity. By inserting the approximation above and neglecting the temperature gradients 
we rewrite the evolution equation for /(r, v,t) as: 


+ V • V - ■ ^)/(r,v,t) - • (v - u(r,t))n(r,t)(()M(r,v,t) = -z/[/(r,v,t) - n(r,t)(/)M(r, v,t)] 


'dt 


kBT 


( 25 ) 












It is interesting to analyze the structure of the self consistent field which appears in the equation of evolution and is 
obtained by using (113 in m- One finds that the field is given by the sum of different forces, dissipative and non 
dissipative: 


$(r, t) = F™^(r, t) + F”"'=(r, t) + F^(r, t). (26) 

where the first term represents the gradient of the potential of mean force, that is the average force of the remaining 
particles on a particle at r, the second term describes the viscous force due to the presence of velocity gradients, 
whereas the last one is due to the existence of thermal gradients. All terms can be identified and computed self- 
consistently once we know the fields n, u, T. Only the first force is non vanishing at equilibrium. For such a reason 
in the DDFT, which only considers the equilibrium free energy, the effective force is purely non dissipative. In fact 
F™/(r,t) may only describe forces which have an equilibrium counterpart entropic, depletion, electrostatic. Van der 
Waals, etc, but not velocity dependent forces or forces due to thermal gradients. 

Using the explicit form of the RET collision dTS]) the following expressions for the different forces are derived: 

{r,t) = —ksTa^ J dkkig 2 {r,r + ak,t)n{r + ak,t). (27) 

where one can show that for small density gradients: F'"^(r,t) = —\7g.exc{^,t). The expression of the viscous force 
reads 


^OTsc(r, ^) _ 2(7^ J dfcfci52(r, r + crfc, t)n{r + ak, t)kj{uj{r + ak, t) — Uj{r, t)) 

and the force due to the presence of thermal gradients is 

F^{r,t) = -^ / dkkig 2 ir,r + ak,t)n{r + ak,t)kB[T{r + ak,t) - T{r,t)]. (28) 

Notice that all these forces are expressed as convolutions involving g 2 (r, r', t), n(r, t), u(r, f), T(r, t) . From the equa¬ 
tions above one can derive the following explicit expressions for the collisional contributions to the bulk dynamical 
viscosity coefficient: 


= —\/TTmkBT a^g2{o')r 
15 


(29) 


and to the heat conductivity: 


y/rrurk^g2{(j)n^ 

3 m 

which correspond to the formulae proposed by Longuet-Higgins and Pople [s^. 

Higher order corrections which better describe the density dependence of transport coefficient have also been 
included by employing an extended approximation which takes into account the dependence of the relaxation time v 
on the hydrodynamic modes [s^. 

Before closing this section we would like to mention the important work performed in a spirit similar to our work 
by L.S. Luo and coworkers (ssl. and more recently by Baskaran and Lowengrub (3^ . 


V. APPLICATIONS OF THE WDLBM TO NEUTRAL AND CHARGED MIXTURES 

In the present section we briefly discuss recent extensions and applications of the self-consistent dynamical method, 
to binary neutral or to ternary charged mixtures comprised of positively and negatively charged hard spheres carrying 
point like charges (the ions) plus neutral spheres representing the solvent [^. Each species, denoted by a, has mass 
m“, diameter cr““. In the case of the binary mixture there are no charges (valence z°‘ = 0) and there are only two 
components, while for the Coulomb case the two ionic components carry charges 2 :“e, e being the proton charge . In 
the latter case the index a = 0, identifies the solvent , whose valence is zero, while a = ± identifies the two oppositely 
charged ionic species interacting through a uniform medium of constant dielectric permittivity. 
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We describe the system by the following set of Enskog-like equations governing the evolution one-particle phase 
space distributions /“(r,v,t) of each species [s^ : 

d F“('r) d 

-r(r,v,t) + V . Vr(r,v,t) + ^ • ^r(r,v,t) = 

- ^ • (v- u(r,t))n“(r,r)(;i“(r,v,t) = — • —r{Y,v,t) (30) 

kbi 'rn°‘ av 

The BGK term is modified with respect to eq (1251) and contains the distributions functions 0“ and ()>“ which have 
the following representations (see ref. for details): 


‘(i',v,i) = [ 


2'KkBT 


13/2 


exp -- 


m“(v — u(r, t)y 


2kBT 


( 31 ) 


and 


m“(u“(r,t) — u(r,t)) • (v — u(r,t)) 


l(r,v,t) = <()“(r,v,t)|l 

to“ ^m“[(u“(r,t)-u(r,t))-(v-u(r,t))]^ ^ 

■^(- k^ - ir,t)-uiY,t)) 


kBT 
1 2 


(32) 


In the case of a one-component fluid there is no difference between <()“ and <()“, since the velocities u“ and u coincide 
and the standard BGK approximation involves the difference between the distribution /“ and the local Maxwellian 
n°‘{Y,t)(j)°‘{Y,v,t). The reason to use the modified distributions (IMll instead of (1311) is to obtain the correct mutual 
diffusion and hydrodynamic properties starting from (I5(I|) . By taking a simple BGK recipe, that is, by setting </>“ =</>“, 
would lead to a double counting of the interactions on the diffusive properties. 

The effective field acting on each species, $“(r, t), stems from interparticle forces and can now be separated in the 
case of isothermal systems into three contributions by generalizing eq. (|M)l (4ll | as 

$“(r, t) = F“’™l(r, t) + F“’‘^™s(r, t) + F“’“"'=(r, t). (33) 

The new term F“’‘^’’“s is a frictional force resulting from the fact that the velocities of each species can be different 
from the barycentric fluid velocity and depends linearly on their relative velocities: 

F“’‘'™9(r, t) = -J2 0(u“(r, t) - u^(r, t)) (34) 

P 


where 7 “-^ is the inhomogeneous friction tensor given by: 


7“^(r,t) =2((T“^)2y^^^Ji? J dk4fcj<?“'^(r,r-Pa“^k,t)n'5(r-Pa“^k,t), 


(35) 


where is the reduced mass , (t“^ the arithmetic average of the diameters for the colliding pair 

and g°‘^ is the pair correlation at contact for species a and /3 whose value is obtained from the extension of the 
Garnahan-Starling equation to mixtures (d^ - ld^ l. We have shown (see ref. 0 ) that when there exists a disparity of 
masses and concentrations between the two species of a binary mixture, the heavier and diluted component can be 
treated as an assembly of colloidal particles moving in a fluid assimilated to a solvent, whose effect is to exert a drag 
force on the colloidal particles. One eliminates the fluid variables relative to the light component from the description 
of the composite system and obtain a closed equation for the heavy particles only . To achieve that, one assumes that 
the dynamical properties of the heavy particles evolve on a time-scale much longer than the characteristic time-scale 
of the light particles. Calling n'^(r,t) the number density of colloidal particles and n®(r,t) the analogous quantity for 
the solvent particles in the dilute limit, jn^ << 1 , one obtains the following advection-diffusion equation: 


dtrf{Y, t) + V • (u(r, t)rf{Y, t)) 



n^(r,t)(v/i'=(r,t) 



( 36 ) 
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where for a bulk system the friction coefficient, 7 , is given by the expression: 


7 8n^ y'Trm^kBTa^^gcsicTcs) 

and the local chemical potential of the colloidal particles is determined by the colloidal-colloidal and by the colloidal- 
solvent interactions, /r^(r, t), (see refs, ): 

VAi^(r,t) = fcBTVlnn^(r,t) - F‘=’™'^(r, t). 


We, now, observe that ea. (|36)l is formally identical to a DDFT equation for the c species in a velocity field u(r,t). 
As —j> 0 the gradient of the chemical potential approaches the ideal gas value, kBT\7n°{r,t), so that eq. (|36ll 
becomes a linear advection-diffusion equation for the field n°, with a diffusion coefficient given by: 


D = 


kBT 

1 


(38) 


to be interpreted as a fluctuation-dissipation relation between 7 and D . 


A. Charged mixtures 

In the case of charged mixture one has to take separately into account the electrostatic potential generated 

by the charge distribution = e(n+(r,t) — (where are the zeroth velocity moments of /^) and by 

the fixed charges located on the pore surfaces and on the electrodes and satisfies the Poisson equation: 

V2?/;(r,t) =(39) 
e 

with boundary conditions —Vifj{r,t) ■ n = S(r)/e at the confining surfaces, where i;(r) is the surface charge density 
sitting on the boundaries and h is the unit vector normal to the surface. 

In the limit of slowly varying fields the transport equations (1301) must reproduce the equations of Electrokinetics. 
In order to recover the equation describing the coupling between diffusion and drift induced by the presence of electric 
field, the so called Poisson-Nernst-Planck equation, one uses the momentum balance equation for the species a and 
neglects inertial terms and finds 


dn^{r,t) 

Ft 


+ V-J±(r,t) =0 


(40) 


with 


jf (r,t) = + n^{r,t)u{r,t) (41) 

where 7 ^ represents the drag coefficient due to the frictional force exerted by the fluid on the particles of type a = ±, 
in reason of their drift velocities and the average barycentric velocity u is given by 


u(r,t) 


X)„m“n°^(r,f)u°(r,f) 

m“n“(r,t) 


(42) 


In dilute solutions the charged components are expected to experience a large friction arising only from the solvent 
while a negligible friction from the oppositely charged species, so that we further approximate and the friction can be 
evaluated in uniform bulk conditions to be 


2TTkBT- 




(43) 


with vP being the bulk density of the solvent and the bulk ion-solvent pair correlation function evaluated at 
contact. In stationary conditions, the following approximated force balance is obtained 

VpF (r, t) - F± (r) -b ez± VV’(r, t) « (r, t). 


(44) 
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Finally, using eq. (1341) we obtain an expression for the ionic currents in terms of the microscopic parameters which 
has the same form as the phenomenological Planck-Nernst current (HD) , with the full chemical potential gradient 
replacing the ideal gas chemical potential gradient, ksTW \nn°‘ The total electric charge density current is 

_^ 

Je = - ^ ^^n^(r,t)V/r=‘=(r,t) + UezE + pe{r,t)u{r,t) (45) 


where the zero frequency electric conductivity aei for a uniform system is given by the Drude-Lorentz-like formula: 


^el — ^ 


T 


T 


(46) 


showing that the conductivity is due to collisions with the solvent and decreases as the solvent becomes denser ( 7 =*= 
is an increasing function of •nP = ‘nP = n~) while increases with the number of charge carriers. 

Finally we obtain the total momentum equation 


dtUj{r,t) + uPr,t)\7iUj{r,t) + ez^n^{r,t)\7 

Pm Pm ^ 

- — ^ n“(r,t)(F;“(r)+F“’™^(r,t)+F;“'"«'=(r,t))=0. (47) 

The total kinetic pressure which can be approximated as: 

knTS^jV, Y. n“(r, t) - , (48) 

a 

with n“. Using the result of ref. we can write 

V n“(r, <)F“’”""(r, t) ~ V^u - + rjP)V{V ■ u). (49) 

(j 

a. 


The non-ideal contribution to the shear viscosity is evaluated in uniform bulk conditions to be 

^ ^ ^ \/27ryi“/5fcBT(cr“^)^5“^n^n^, 


(50) 


In conclusion, (lT7l) can be cast in the Navier-Stokes form 

1 - 

dtUj + UiViUj = - S/iPSij — —Vjip + —S/iViUj —VjViUi (51) 

Pm Pm Pm Pm 

with r] = and = ViPid + Sa ^“(i’U)Vi/r“ 3 ,c(r, t). Eqs. (|^ together with the continuity equation 

for each species, equation m and the Poisson equation are sufficient to understand the behavior of an electrolyte 
solution. 

In refs. [H, we have applied these equations to the study of electrokinetics flow in channels of nanometric 
section and obtained the current-voltage relations and the conductance. In addition we have also considered the 
electro-osmosis phenomenon by which an electric field can be used to induce a mass flow in an electrolyte solution 
even in the absence of pressure gradients. 


VI. NUMERICAL SOLUTION OF THE KINETIC EQUATION VIA LBM 

In this section we briefly illustrate how the LBM is applied as a numerical solver to the present kinetic model, which 
we shall discuss by considering only the one component case, while the multicomponent case can be easily deduced. 
Other popular schemes such as BGK or Shan-Chen collisional forms require a minor numerical effort, but on the other 
hand fail to predict non trivial transport coefficients [lTj |. 

The WDLBM consists in integrating directly the kinetic equation equation (1251) for /(r, v,t) using a discretization 
procedure introduced some 20 years ago in fluid dynamics The LBM differs from other methods that are 
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based on the solution of the set of equations for the density and the velocity of the fluid or from the Direct Monte Carlo 
simulation (5lj |. The numerical method we have adopted is a substantial modification of the conventional method 
used in fluid dynamics applications to the presence of hard sphere collisions. 

The strategy behind the LBM is to reduce the number of possible particle spatial positions and microscopic velocities 
V from a continuum to just a finite number of values, Cp, and similarly discretize time into distinct steps. The typical 
evolution equation is a time-explicit integration that reads 

|[+vV/ = D'(/,M) (52) 

where the kernel D' contains both the collisional term, the BGK term and the external force term F • ^/, while M 
represents a generic moment of /. 

We begin by projecting the phase space distribution function over an orthonormal basis spanned by the tensorial 
Hermite polynomials {H^}: 


LXJ ^ 

/(r, V, t) = w(v) ^ (53) 

1=0 

where a;(v) = (27ru|n)“^/^e“'' and, using the orthonormal relation 

J (v)dv = (vT)''~'''^Sim.Sg£ (54) 

the moments Ma'^ (r, t) can be obtained by: 

J /(r,v,t)iJ«(v)dv (55) 

The exact infinite series representation given by eq. (1531) is approximated by a function /(r, v, t) obtained by retaining 
in eg. (1531) only terms up to Z = K. Using the Gauss-Hermite quadrature formula in order to evaluate the expansion 
coefficients, and using the nodes, Cp, and the weights Wp of a quadrature of order 2G > K one obtains 

G 

= ^/p(r,t)ijW(cp) (56) 

p—0 

with 

/p(r, t) = /(r, Cp, <)-^ (57) 

w(Cp) 

and /(r,v,t) —>• fp{r,t) is the truncation of order K of the expansion of eq. (l53l) . 

The distribution function is replaced by an array of Q populations (19 in the standard three dimensional case) and 
the distribution function /(r, v, t) —?► /p(r, t)- The propagation of the populations is achieved via a time discretization 
to first order and a forward Euler update: 


dtfpir, t)+Cp- drfpir, t) ~ 


fp(v -k CpAt, t + At)- /p(r, t) 
At 


(58) 


where At is the time-step. 

By the same token, we consider the expansion of the collisional kernel. 


D(r, V, t) = a;(v) Y (’*’ ('"^ 


(59) 


1=0 

The operational version of the collisional term is provided by 


Dp(r,t) = D(r,Cp,t) 


w(cp) 


(60) 
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From these transformations, the evolution equation of the new representation is given by the following updating 
scheme 


/p(r + Cp At, t + At) = /p(r, t) + flp(r, t) At (61) 

where 

ftp(r,t) = :/(/p«(r,t) -/p(r,t)) + 5'p(r,t) (62) 

The explicit form of the r.h.s. of eqs. (IMli-dMl) reads 

= Wp n(r,t)+n(r,t)u(r,t) •'H^^)(cp)+n(r,t)u(r,t)u(r,t) :'H^^)(cp) (63) 

5p(r,t) = u;p{n(r,t)(F-/(r,t)+F™-(r,t)+F“‘(r)). [?tW(cp)+2?t(2)(cp).u(r,t)]} (64) 

where 'H^^^(cp) = ^ and 'H^^^(cp) = being a vector and a tensor of rank two, respectively, and I is the unit 

tensor. 

Once the populations fp are known, they are used to compute hydrodynamic moments, entering both the equilibrium 
and in sampling the macroscopic evolution. The fluid density, momentum current and local temperature read 

n{r,t) = ^/p(r,t) 

P 

n(r,t)u(r,t) = ^/p(r,t)cp 

P 

;(r,t)r(r,t) = ^/p(r,t)c2 


3fc_B 


m 


(65) 


The dynamics of the charged multicomponent system can be solved numerically in the same framework by general¬ 
izing the procedure described above to three species. The novelty consists in the presence of the electric field t/'(r, t) 
which must be determined from the charge distribution by solving the Poisson equation for the electrostatic poten¬ 
tial generated by the mobile and surface charges. The determination of '(/’(r) can be achieved by using a successive 
over-relaxation method, while the speed of convergence is greatly enhanced by employing a Gauss-Seidel checker¬ 
board scheme in conjunction with Chebychev acceleration j52j[ . Neumann boundary conditions on the gradient of the 
electrostatic potential are imposed at the wall surface 

fl • VV'IrGS = (66) 

e 

where n is the normal to the surface, S. 


VII. CONCLUSIONS 

In summary, we reviewed recent advances in kinetic theory and modeling applied to the transport of molecular 
liquids confined in very small spaces. We formulated the model and then considered the series of approximations 
needed in order to achieve a workable numerical scheme to be used under generic confinement conditions. 

At first, we split the collision operator into hydrodynamic and non-hydrodynamic (i.e. purely kinetic) contributions. 
The non-hydrodynamic contribution is handled by the BGK ansatz, which gives rise to the ideal gas thermodynamics 
and to non-collisional stress terms. The non hydrodynamic part instead is treated by performing the integrals featuring 
in the collision operator by using a parametric form of the phase space distribution /. Since the collision operator 
involves convolutions with the configurational pair correlation function we used the so-called adiabatic approximation 
to determine it. The resulting equation involves only the single particle distribution function and its first few moments, 
and can be solved by via the Lattice Boltzmann method. Such strategy gives access to the complete hydrodynamic 
information without the need to solve the hydrodynamic equations separately. 
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